#!/usr/bin/perl -w
# given a seed sequence, a parameter file, a codon frequency table, and a set
# of this restriction enzyme sites (forbidden subsequences), this script 
# generates diverse protein-coding sequences.
use strict;
use warnings;

if (@ARGV != 0) {
  die "USAGE: $0\n";
}

my $seed = getSingleFastaSequence("protein.fasta");
my $effector = getSingleFastaSequence("effector.fasta");
my @parameters = getParameters("params.txt");
my %codonFreqs = getCodonFreqs("codonfreq.txt");
my @sites = getRestrictionSites("sites.fasta");
my $quantity = $parameters[0];
my $outputfile = "output.fasta";
open OUTPUT, ">$outputfile" or die "couldn't open $outputfile: $!";

for (my $i=0; $i<$quantity; $i++) {
  my $sequence = getCodingSequence($seed, \@parameters, \%codonFreqs, \@sites, $effector);
  for my $rs (@sites) {
    if($sequence =~ m/$rs/i) {
      print "Sequence matches restriction site. Regenerating...\n";
      $i--;
      next;
    }
  }
  print "Sequence ", $i+1, "...\n";
  if ($sequence eq 'noseq') {
    print "Failed. Re-Attempting.\n";
    $i--;
  } else {
    print OUTPUT ">generated sequence ", $i+1, "\n";
    $sequence =~ s/(.{1,80})/$1\n/g;
    print OUTPUT $sequence;
  }
}

sub getCodingSequence {
  my ($seed, $parameters, $codonFreqs, $sites, $effector) = @_;
  $seed = uc $seed;
  my @parameters = @{$parameters};
  my %codonFreqs = %{$codonFreqs};
  my @sites = @{$sites};
  # codon to amino acid mapping; 'X' signifies stop codon
  my %codonTable = (
   "AAA" => "K", "AAC" => "N", "AAG" => "K", "AAU" => "N",
   "ACA" => "T", "ACC" => "T", "ACG" => "T", "ACU" => "T",
   "AGA" => "R", "AGC" => "S", "AGG" => "R", "AGU" => "S",
   "AUA" => "I", "AUC" => "I", "AUG" => "M", "AUU" => "I",

   "CAA" => "Q", "CAC" => "H", "CAG" => "Q", "CAU" => "H",
   "CCA" => "P", "CCC" => "P", "CCG" => "P", "CCU" => "P",
   "CGA" => "R", "CGC" => "R", "CGG" => "R", "CGU" => "R",
   "CUA" => "L", "CUC" => "L", "CUG" => "L", "CUU" => "L",

   "GAA" => "E", "GAC" => "D", "GAG" => "E", "GAU" => "D",
   "GCA" => "A", "GCC" => "A", "GCG" => "A", "GCU" => "A",
   "GGA" => "G", "GGC" => "G", "GGG" => "G", "GGU" => "G",
   "GUA" => "V", "GUC" => "V", "GUG" => "V", "GUU" => "V",

   "UAA" => "X", "UAC" => "Y", "UAG" => "X", "UAU" => "Y",
   "UCA" => "S", "UCC" => "S", "UCG" => "S", "UCU" => "S",
   "UGA" => "X", "UGC" => "C", "UGG" => "W", "UGU" => "C",
   "UUA" => "L", "UUC" => "F", "UUG" => "L", "UUU" => "F"
  );
  # 0 to 100 percent ID divided by sequence length
  my @expectedFractionGapsPerSeqLen = (0.100600600600601, 0.103103103103103, 0.103603603603604, 0.0991991991991992, 0.1003003003003, 0.0968968968968969, 0.100600600600601, 0.0978978978978979, 0.1003003003003, 0.0924924924924925, 0.0947947947947948, 0.08998998998999, 0.0912912912912913, 0.0876876876876877, 0.0871871871871872, 0.0844844844844845, 0.0771771771771772, 0.0750750750750751, 0.0695695695695696, 0.0696696696696697, 0.0692692692692693, 0.0614614614614615, 0.0611611611611612, 0.062962962962963, 0.0531531531531532, 0.0541541541541542, 0.0517517517517518, 0.0452452452452452, 0.0467467467467468, 0.044044044044044, 0.0404404404404404, 0.04004004004004, 0.0352352352352352, 0.0397397397397397, 0.0351351351351351, 0.0309309309309309, 0.0334334334334334, 0.0294294294294294, 0.027027027027027, 0.0287287287287287, 0.0232232232232232, 0.0242242242242242, 0.0232232232232232, 0.0196196196196196, 0.0236236236236236, 0.0182182182182182, 0.018018018018018, 0.0171171171171171, 0.0154154154154154, 0.0168168168168168, 0.0156156156156156, 0.0125125125125125, 0.0133133133133133, 0.0116116116116116, 0.0106106106106106, 0.0112112112112112, 0.00940940940940941, 0.00940940940940941, 0.00850850850850851, 0.00940940940940941, 0.00790790790790791, 0.00730730730730731, 0.00690690690690691, 0.00720720720720721, 0.0058058058058058, 0.00600600600600601, 0.00560560560560561, 0.00530530530530531, 0.0038038038038038, 0.0039039039039039, 0.0039039039039039, 0.0051051051051051, 0.0024024024024024, 0.00390390390390391, 0.0023023023023023, 0.0021021021021021, 0.0019019019019019, 0.0019019019019019, 0.0028028028028028, 0.0016016016016016, 0.0017017017017017, 0.0012012012012012, 0.0013013013013013, 0.0014014014014014, 0.0015015015015015, 0.000900900900900901, 0.0005005005005005, 0.001001001001001, 0.000700700700700701, 0.0005005005005005, 0.0005005005005005, 0.0003003003003003, 0.0002002002002002, 0.0001001001001001, 0.0001001001001001, 0, 0.0001001001001001, 0, 0, 0, 0);

  my @forbidden = @sites;
  push @forbidden, $effector;
  for (my $i=0;$i<@forbidden;$i++) {
    $forbidden[$i] = uc $forbidden[$i];
  }
  my %reverseCodonTable = getReverseCodonTable(%codonTable);
  my $qty=$parameters[0];
  my $fid=$parameters[1];
  my $gapPercent = $parameters[2];
  my $gapFraction = $gapPercent/100;
  my $flooredFID = int($fid*100)/100;
  # if the maximum percent of allowed gaps is greater than the percent we expect, the percent of allowed gaps must decrease by increasing the fraction id
  while (($flooredFID <= 1) and ($gapFraction < $expectedFractionGapsPerSeqLen[$flooredFID*100])) {
    $flooredFID+=0.01;
    print "Had to raise FID. New Value=$flooredFID\n";
  }
  if (int($fid*100)/100 != $flooredFID) {
    $fid=$flooredFID;
  }
  my %nonStopCodonProbs = getNonStopCodonProbs(%codonFreqs);
  my $longestForbiddenLength = getLongest(@forbidden);
  my $initSeq = 
   getInitialSequence($seed,\%reverseCodonTable,\%codonFreqs,\@forbidden,$longestForbiddenLength);
  my $mutatedSeq = 
   mutate($initSeq,\%codonTable,\%nonStopCodonProbs,\@forbidden,$fid,\%codonFreqs, $longestForbiddenLength);
  my $formattedSeq = 
   formatSequence($mutatedSeq,\@forbidden,\%codonFreqs, $longestForbiddenLength);
  my $translatableSeq = $formattedSeq;
  # get the protein sequence
  $translatableSeq =~ tr/acgt/ACGT/;
  $translatableSeq =~ tr/T/U/;
  my @codonSeq = ($translatableSeq =~ m/.{3}/g);
  pop(@codonSeq);
  my @aminoSeq = ();
  for(my $i=0; $i<@codonSeq; $i++) {
    $aminoSeq[$i] = $codonTable{$codonSeq[$i]};
  }
  my $proteinSeq = join("",@aminoSeq);
  my @alignment = needlemanWunsch($seed, $proteinSeq);
  # print $alignment[1], "\n";
  $alignment[1] =~ s/[^-]//g;
  my $observedGapFraction = length($alignment[1])/length($seed);
  print "observed gap fraction:$observedGapFraction\n";
  if ($observedGapFraction > $gapFraction) {
    return "noseq";
  } else {
    return $formattedSeq;
  }
}

sub getLongest {
  my (@inputs) = @_;
  my $longest = length($inputs[0]);
  foreach my $input (@inputs) {
    if (length($input) > $longest) {
      $longest = length($input);  
    }
  }
  return $longest;
}

# reverse the codon table so that amino acids are the keys and lists of codons
# are the values
sub getReverseCodonTable {
  my (%codonTable) = @_;
  my %reverseCodonTable;
  foreach my $codon (keys %codonTable) {
    my $amino = $codonTable{$codon};
    if (not defined($reverseCodonTable{$amino})) {
      $reverseCodonTable{$amino} = [$codon];
    } else {
      push @{$reverseCodonTable{$amino}}, $codon;
    }
  }
  return %reverseCodonTable;
}

# get protein.fasta
sub getSingleFastaSequence {
  my ($fastafile) = @_;
  # seed protein sequence
  open INPUT, $fastafile or die "couldn't find $fastafile: $!";

  my $seqLabel;
  my $seq = "";
  while(my $line = <INPUT>) {
    chomp $line;
    if($line =~ m/^\s*$/) {
      next;
    } elsif ($line =~ m/^>/) {
      $line =~ s/^>//;
      $seqLabel = $line;
    } else {
      $line =~ s/\s//g;
      $seq .= $line;
    }
  }

  close INPUT;
  return $seq;
}
# get the number of sequences, target percent ID between seed and generated sequences,
# and maximum percentage of allowed gaps when generated sequences are aligned 
# to seed sequence
sub getParameters {
  my ($paramfile) = @_;
  open INPUT, $paramfile or die "couldn't find $paramfile: $!";

  my @params;
  while(my $line = <INPUT>) {
    chomp $line;
    if ($line =~ m/^\s*$/) {
      next;
    } else { 
      $line =~ s/^.*?([-e\.\d]+).*$/$1/;
      push @params, $line;
    }
  }
  close INPUT;

  # turn fraction ID into percent ID:
  $params[1] /= 100;
  return @params;
}

# Codon frequencies
sub getCodonFreqs {
  my ($codonfreqfile) = @_;
  open INPUT, $codonfreqfile or die "couldn't find $codonfreqfile: $!";

  my %codonFreqs;
  while(my $line = <INPUT>) {
    if($line =~ m/^\s*$/) {
      next;
    } else {
      $line =~ s/^[ \t]*//;
      my @pair = split(/\s+/, $line);
      $codonFreqs{$pair[0]}=$pair[1];
    }
  }
  close INPUT;
  return %codonFreqs;
}

# get the probabilities for each codon that isn't a stop codon
sub getNonStopCodonProbs {
  my (%codonFreqs) = @_;
  my $total=0;
  my @codons = (keys %codonFreqs);
  my %nonStopCodonProbs;
  foreach my $codon (@codons) {
    if( $codon !~ m/UGA|UAA|UAG/i) {
      $total+=$codonFreqs{$codon};
    }
  }
  foreach my $codon (@codons) {
    if( $codon !~ m/UGA|UAA|UAG/i) {
      $nonStopCodonProbs{$codon} = $codonFreqs{$codon}/$total;
    }
  }
  return %nonStopCodonProbs;
}

# get the restriction enzyme sites
sub getRestrictionSites {
  my ($sitesfile) = @_;
  open INPUT, $sitesfile or die "couldn't find $sitesfile: $!";

  my $seqCount=0;
  my @siteLabels;
  my @sites;
  while(my $line = <INPUT>) {
    chomp $line;
    if($line =~ m/^\s*$/) {
      next;
    } elsif($line =~ m/^>/) {
      $line =~ s/^>//;
      $seqCount++;
      push @siteLabels, $line;
      $sites[$seqCount-1]="";
    } else {
      $line =~ s/\s//g;
      $sites[$seqCount-1] .= $line;
    }
  }
  close INPUT;

  my @forbidden;
  push @forbidden, (@sites);
  foreach my $site (@sites) {
    push @forbidden, getReverseComplement($site);
  }

  return @forbidden;
}

# generate the reverse complement of a single sequence
sub getReverseComplement { 
  my ($sequence) = @_;
  my %WCPairs = (
      "A" => "T",
      "T" => "A",
      "G" => "C",
      "C" => "G"
  );
  $sequence = uc $sequence;
  $sequence =~ tr/U/T/;
  my @forward = ($sequence =~ m/./g);
  my @rev = reverse @forward;
  my @revCompl;
  for(my $i=0; $i<@rev; $i++) {
    my $letter = $rev[$i];
    my $flip = $WCPairs{$letter};
    push @revCompl, $flip;
  } 
  return join("", @revCompl);
}

# get an initial sequence
# generating the initial sequences linearly makes it faster and easier to
# assert that there are no restriction sites
sub getInitialSequence {
  my ($seed, $reverseCodonTable, $codonFreqs, $sites, $restSiteLen)  = @_;
  my %reverseCodonTable = %{$reverseCodonTable};
  my %codonFreqs = %{$codonFreqs};
  my @codons = (keys %codonFreqs);
  my $sequence = "";
  for(my $j=0; $j<length($seed); $j++) {
    my $amino = uc substr $seed, $j, 1;
    my @codons = @{$reverseCodonTable{$amino}};
    my %freqs;
    foreach my $cod (@codons) {
      $freqs{$cod}=$codonFreqs{$cod};
    }
    my $trialSeq = $sequence . randcod(%freqs);
    if(containsForbiddenSequence($trialSeq, $sites, $restSiteLen, $j*3) eq 'true') {
      # delete the previous 3 codons and restart
      $sequence=substr $sequence, 0, (($j-2>0)?$j-2:0)*3;
      $j=($j-3 > -1) ? $j-3 : -1;
      # print STDERR "Produced restriction sequence, regenerating a few previous codons...\n";
    } else {
      $sequence = $trialSeq;
    }
  }
  return $sequence;
}

sub mutate {
  my ($initSeq, $codonTable, $nonStopCodonProbs, $sites, $fid, $codonFreqs,$longestRestLen) = @_;
  my $seedLength = length($initSeq)/3;
  my $sequence=$initSeq;
  my $qtyMutatedSites=int((1-$fid)*($seedLength-1));
  my %counts;
  foreach my $codon (keys %{$nonStopCodonProbs}) {
    $counts{$codon}=0;
  }
  foreach my $codon ($initSeq =~ m/.{3}/g) {
    $counts{$codon}++;
  }
  # Strategy: 
  my (@bagOfCodonSites) = 1 .. ($seedLength-1);
  # for each necessary mutation
  for (my $j=0; $j<$qtyMutatedSites; $j++) {
    my $failure;
    do {
      # 1. select a random site in the sequence to insert the codon at.
      my $draw = splice @bagOfCodonSites, int(rand(@bagOfCodonSites)), 1;
    #   print STDERR "Remaining sites: ", join(" ", @bagOfCodonSites), "\n";
    #   print STDERR "Selected site: $draw\n";
      # 2. generate a random codon according to a distribution that will bring the sequence closer to the background distribution.
      my $prevCodon = substr $sequence, $draw*3, 3;
    #  print STDERR "Previous codon: $prevCodon\n";
      my $codon = 
       randcod(idealDistrib($nonStopCodonProbs, $codonFreqs, \%counts, $seedLength, $prevCodon, $codonTable));
    #   print STDERR "New codon: $codon\n";
      # 3. assert that the sequence does not match any forbidden sequence.
      my $trialSeq = $sequence;
      substr ($trialSeq, $draw*3, 3, $codon);
    #   print STDERR "Original sequence: ", join(" ",($sequence=~m/.{3}/g)), "\n";
    #   print STDERR "Trial sequence:    $trialSeq\n";
      # 4. if it does, just put the position back in the bag and pick another one
      $failure = containsForbiddenSequence($trialSeq, $sites, $longestRestLen, $draw*3);
      if($failure eq 'true') {
     #   print STDERR "Produced restriction sequence; choosing different site...\n";
        push @bagOfCodonSites, $draw;
      } else {
        $sequence = $trialSeq;
        $counts{$prevCodon}--;
        $counts{$codon}++;
      }
    } while ($failure eq 'true');
  }
  return $sequence;
}

sub formatSequence {
  my ($initSeq, $sites, $codonFreqs, $longestRestLen) = @_;
  my $sequence = $initSeq;
  # format start codon
  $sequence = lc $sequence;
  my $startCodon = substr $sequence, 0, 3;
  $startCodon = uc $startCodon;
  substr $sequence, 0, 3, $startCodon;
  # get stop codon
  my $failure;
  do {
    my $stopCodon = generateStopCodon($codonFreqs);
    $stopCodon = uc $stopCodon;
    my $trialSeq = $sequence . $stopCodon;
    $failure = containsForbiddenSequence($trialSeq, $sites, $longestRestLen,length($sequence));
    if ($failure eq 'false') {
      $sequence = $trialSeq;
    }
  } while( $failure eq 'true' );
  $sequence =~ tr/uU/tT/;
  return $sequence;
}

# return a hash containing the probability distribution that will be used
# to generate the next codon
sub idealDistrib {
  my ($nonStopCodonProbs,$codonFreqs,$counts,$qty,$prevCodon,$codonTable) = @_;
  my %nonStopCodonProbs = %{$nonStopCodonProbs};
  my %codonFreqs = %{$codonFreqs};
  my %counts = %{$counts};
  my %codonTable = %{$codonTable};
  my %neededCodons;
  foreach my $codon (keys %nonStopCodonProbs) {
    my $numCodons = ($nonStopCodonProbs{$codon}*$qty-$counts{$codon})/(1-$nonStopCodonProbs{$codon});
    if ($numCodons > 0 and $codonTable{$prevCodon} ne $codonTable{$codon}) {
      $neededCodons{$codon} = $numCodons;
    }
  }

  return %neededCodons;
}

sub containsForbiddenSequence {
  my ($inseq,$forbidden,$longestLen,$firstBasePos) = @_;
  $inseq =~ tr/uU/tT/;
  my @forbidden = @{$forbidden};
  my $lower = $firstBasePos-$longestLen;
  $lower = ($lower > 0 ? $lower : 0);
  my $upper = $firstBasePos + 3 + $longestLen;
  $upper = ($upper < length($inseq)-1 ? $upper : length($inseq)-1);
  my $subSeq = substr $inseq, $lower, $upper-$lower+1;
  foreach my $rs (@forbidden) {
    $rs =~ tr/uU/tT/;
    if ( $inseq =~ m/$rs/i ) {
      return 'true';
    }
  }
  return 'false';
}

# generate a random codon according to frequency distribution %codonFreq
# probabilities should be changed to reflect initial sequence's bias
sub randcod {
  my (%comp) = @_;
  my $sum = 0;
  foreach my $val (values %comp) {
    $sum += $val;
  }
  my ($rnum, $total) = rand($sum), 0;
  foreach my $cod (keys %comp) {
    $total += $comp{$cod};
    if ($rnum < $total) {
      return $cod;
    }
  }
}

sub generateStopCodon {
  my ($codonFreqs) = @_;
  my %codonFreqs = %{$codonFreqs};
  my %stopCodonFreqs;
  foreach my $codon (keys %codonFreqs) {
    if($codon =~ m/UGA|UAG|UAA/i) {
      $stopCodonFreqs{$codon} = $codonFreqs{$codon};
    }
  }
  return randcod(%stopCodonFreqs);
}

=exanple
my @align = needlemanWunsch("adamwest","admwst");

for (my $i=0; $i<=length($align[0]); $i+=80){
  print substr($align[0],$i,80)."\n";
  print substr($align[1],$i,80)."\n";
  print "\n";
}
=cut

# the needleman-wunsch pairwise alignment algorithm, courtesy of the be131 homework solutions
sub needlemanWunsch {
  my (@seqs) = @_;
  my @s = ();
  my @tb = ();
  my $tempMax = 0;
  my $pair = '';
  my $nextI = 0;
  my $nextJ = 0;
  my $nextPos = '';
  my $seq1 = '';
  my $seq2 = '';
  my $flagEnd =0;

  for(my $i=0;$i<2;$i++) {
    $seqs[$i] = uc $seqs[$i];
  }

  my $gapPenalty = -4;
  my @matrixChars = ("A", "R", "N", "D", "C", "Q", "E", "G", "H", "I", "L", "K", "M", "F", "P", "S", "T", "W", "Y", "V", "B", "Z", "X");
  my %preMatrix = 
    ( "A" => [4, -1, -2, -2, 0, -1, -1, 0, -2, -1, -1, -1, -1, -2, -1, 1, 0, -3, -2, 0, -2, -1, 0],
      "R" => [-1, 5, 0, -2, -3, 1, 0, -2, 0, -3, -2, 2, -1, -3, -2, -1, -1, -3, -2, -3, -1, 0, -1],
      "N" => [-2, 0, 6, 1, -3, 0, 0, 0, 1, -3, -3, 0, -2, -3, -2, 1, 0, -4, -2, -3, 3, 0, -1],
      "D" => [-2, -2, 1, 6, -3, 0, 2, -1, -1, -3, -4, -1, -3, -3, -1, 0, -1, -4, -3, -3, 4, 1, -1],
      "C" => [0, -3, -3, -3, 9, -3, -4, -3, -3, -1, -1, -3, -1, -2, -3, -1, -1, -2, -2, -1, -3, -3, -2],
      "Q" => [-1, 1, 0, 0, -3, 5, 2, -2, 0, -3, -2, 1, 0, -3, -1, 0, -1, -2, -1, -2, 0, 3, -1],
      "E" => [-1, 0, 0, 2, -4, 2, 5, -2, 0, -3, -3, 1, -2, -3, -1, 0, -1, -3, -2, -2, 1, 4, -1],
      "G" => [0, -2, 0, -1, -3, -2, -2, 6, -2, -4, -4, -2, -3, -3, -2, 0, -2, -2, -3, -3, -1, -2, -1],
      "H" => [-2, 0, 1, -1, -3, 0, 0, -2, 8, -3, -3, -1, -2, -1, -2, -1, -2, -2, 2, -3, 0, 0, -1],
      "I" => [-1, -3, -3, -3, -1, -3, -3, -4, -3, 4, 2, -3, 1, 0, -3, -2, -1, -3, -1, 3, -3, -3, -1],
      "L" => [-1, -2, -3, -4, -1, -2, -3, -4, -3, 2, 4, -2, 2, 0, -3, -2, -1, -2, -1, 1, -4, -3, -1],
      "K" => [-1, 2, 0, -1, -3, 1, 1, -2, -1, -3, -2, 5, -1, -3, -1, 0, -1, -3, -2, -2, 0, 1, -1],
      "M" => [-1, -1, -2, -3, -1, 0, -2, -3, -2, 1, 2, -1, 5, 0, -2, -1, -1, -1, -1, 1, -3, -1, -1],
      "F" => [-2, -3, -3, -3, -2, -3, -3, -3, -1, 0, 0, -3, 0, 6, -4, -2, -2, 1, 3, -1, -3, -3, -1],
      "P" => [-1, -2, -2, -1, -3, -1, -1, -2, -2, -3, -3, -1, -2, -4, 7, -1, -1, -4, -3, -2, -2, -1, -2],
      "S" => [1, -1, 1, 0, -1, 0, 0, 0, -1, -2, -2, 0, -1, -2, -1, 4, 1, -3, -2, -2, 0, 0, 0],
      "T" => [0, -1, 0, -1, -1, -1, -1, -2, -2, -1, -1, -1, -1, -2, -1, 1, 5, -2, -2, 0, -1, -1, 0],
      "W" => [-3, -3, -4, -4, -2, -2, -3, -2, -2, -3, -2, -3, -1, 1, -4, -3, -2, 11, 2, -3, -4, -3, -2],
      "Y" => [-2, -2, -2, -3, -2, -1, -2, -3, 2, -1, -1, -2, -1, 3, -3, -2, -2, 2, 7, -1, -3, -2, -1],
      "V" => [0, -3, -3, -3, -1, -2, -2, -3, -3, 3, 1, -2, 1, -1, -2, -2, 0, -3, -1, 4, -3, -2, -1],
      "B" => [-2, -1, 3, 4, -3, 0, 1, -1, 0, -3, -4, 0, -3, -3, -2, 0, -1, -4, -3, -3, 4, 1, -1],
      "Z" => [-1, 0, 0, 1, -3, 3, 4, -2, 0, -3, -3, 1, -1, -3, -1, 0, -1, -3, -2, -2, 1, 4, -1],
      "X" => [0, -1, -1, -1, -2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -2, 0, 0, -2, -1, -1, -1, -1, -1]);

  my %simMatrix;
  for my $amino (keys %preMatrix) {
    for(my $i=0;$i<@matrixChars;$i++) {
      $simMatrix{$amino . $matrixChars[$i]} = $preMatrix{$amino}[$i];
    }
  }

#THIS SECTION PERFORMS THE NEEDLEMAN-WUNSCH ALGORITHM

#initialize the array with gap penalties
  $s[0][0] = 0;
  $tb[0][0] = 'END';
  for (my $i=1; $i<=length($seqs[1]); $i++){
    $s[$i][0] = $i * $gapPenalty;
    $tb[$i][0] = 'U';
  }
  for (my $j=1; $j<=length($seqs[0]); $j++){
    $s[0][$j] = $j * $gapPenalty;
    $tb[0][$j] = 'L';
  }

#perform alignment algorithm, store scores in array @s and traceback
#info in array @tb
  for (my $i=1; $i<=length($seqs[1]); $i++){
    for (my $j=1; $j<=length($seqs[0]); $j++){
      $pair = substr($seqs[1],$i-1,1) . substr($seqs[0],$j-1,1);
      $s[$i][$j] = $s[$i-1][$j-1] + $simMatrix{$pair};
      $tb[$i][$j] = 'D';
      $tempMax = $s[$i-1][$j] + $gapPenalty;
      if ($tempMax > $s[$i][$j]) {
        $s[$i][$j] = $tempMax;
        $tb[$i][$j] = 'U';
      }
      $tempMax = $s[$i][$j-1] + $gapPenalty;
      if ($tempMax > $s[$i][$j]) {
        $s[$i][$j] = $tempMax;
        $tb[$i][$j] = 'L';
      }
    }
  }

#perform traceback and print alignment
  $nextI = length($seqs[1]);
  $nextJ = length($seqs[0]);
  $nextPos = $tb[$nextI][$nextJ];
  while ($flagEnd == 0) {
    if ($nextPos eq 'D') {
      $seq1 = substr($seqs[0],$nextJ-1,1) . $seq1;
      $seq2 = substr($seqs[1],$nextI-1,1) . $seq2;
      $nextI--;
      $nextJ--;
    }
    elsif ($nextPos eq 'U') {
      $seq1 = "-" . $seq1;
      $seq2 = substr($seqs[1],$nextI-1,1) . $seq2;
      $nextI--;
    }
    elsif ($nextPos eq 'L') {
      $seq1 = substr($seqs[0],$nextJ-1,1) . $seq1;
      $seq2 = "-" . $seq2;
      $nextJ--;
    }
    elsif ($nextPos eq 'END') {

      $flagEnd = 1;
    }
    $nextPos = $tb[$nextI][$nextJ];
  }

  return ($seq1, $seq2);
}
